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We present a calculation of the scalar field self-force (SSF) acting on a scalar-charge particle in 
a strong-field orbit around a Kerr black hole. Our calculation specializes to circular and equato- 
rial geodesic orbits. The analysis is an implementation of the standard mode-sum regularization 
scheme: We first calculate the multipole modes of the scalar-field perturbation using numerical 
integration in the frequency domain, and then apply a certain regularization procedure to each of 
the modes. The dissipative piece of the SSF is found to be consistent with the fiux of energy and 
angular momentum carried by the scalar waves through the event horizon and out to infinity. The 
conservative (radial) component of the SSF is calculated here for the first time. When the motion is 
^-^ ] retrograde this component is found to be repulsive (outward pointing, as in the Schwarzschild case) 

I , for any spin parameter a and (Boyer-Lindquist) orbital radius tq. However, for prograde orbits we 

' find that the radial SSF becomes attractive (inward pointing) for ro > r'c(a), where Tc is a critical 

I a-dependent radius at which the radial SSF vanishes. The dominant conservative effect of the SSF 

, ^ , in Schwarzschild spacetime is known to be of 3rd post-Newtonian (FN) order (with a logarithmic 

^ • running). Our numerical results suggest that the leading-order FN correction due to the black 

' hole's spin arises from spin-orbit coupling at 3PN, which dominates the overall SSF effect at large 

, ro. In FN language, the change-of-sign of the radial SSF is attributed to an interplay between the 

y—i ■ spin-orbit term (oc —ar^''"^) and the "Schwarzschild" term (oc r^J'^logro). 

(N ; 

I. INTRODUCTION 

o^: 

5^ ' The gravitational two-body problem is extremely difficult to tackle in a general-relativistic context, due to the 
bJO, intrinsic nonlinearities of the theory. However, when one of the two components is much more massive than the other 

the problem simplifies and can sometimes be attacked via black hole perturbation theory. Nature provides us with 

such extreme mass-ratio systems in the form of compact objects inspiraling into massive black holes in galactic nuclei. 
^ Such systems are key targets for the planned space-based gravitational wave detector LISA (Laser Interferometer 

Space Antenna) [l| . Detection of the gravitational waves and accurate extraction of the physical parameters requires 
' precise theoretical templates of the waveforms, which, in turn, necessitate knowledge of the radiative evolution of the 

system. 

" I , The underlying theoretical problem, in its most fundamental form, is that of a pointlike particle orbiting a black 
hole of a much larger mass. The interaction of the particle with its own gravitational field gives rise to a gravitational 
, self- force (GSF), which is responsible in particular for the radiative inspiral. How to calculate this GSF has been 
■ the subject of extensive study over the last decade Q. The fundamental formalism for calculations of the GSF in 
l ' curved spacetime was first laid down by Mino, Sasaki, and Tanaka ^] and independently by Quinn and Wald Q, 
I> . with important later supplements by Detweiler and Whiting Q, Gralla and Wald Q, Pound Q and Harte [1| (See 
Poisson for a review f9|). The resulting equations of motion are known as the MiSaTaQuWa equations. The analogous 
, self-forced equation of motion for the electromagnetic case was derived by DeWitt and Brehme long ago 10] (with 
d ' corrections by Hobbs [ll|) and reproduced more recently using other methods in • Quinn obtained the equivalent 

results for the scalar field self- force (SSF) [Tsj . 

The MiSaTaQuWa equations of motion are hard to implement directly and so they were later recast into forms more 
amenable to practical calculation. One of the standard methods is the mode-sum scheme first introduced in Ref. [l^ . 
Using this method, self force calculations have been performed for a range of problems. These include calculations 
of the SSF for radial infall [l^, circular fl6|. [l7| and eccentric [l^ orbits; the electromagnetic self- force for eccentric 
orbits [l^, and the GSF for radial infall [20], circular [2l|,[2^, and eccentric orbits [23]. More recently, researchers 
have been exploring alternative calculation methods which are based on direct regularization of the self interaction in 
2+1 and 3+1 dimensions (23 - [26| . Common to all calculations presented so far is the fact that they specialize to the 
simpler (but less astrophysically relevant) case where the central object is a non-rotating, Schwarzschild black hole. 

In this paper we open a new front in self force calculations by considering extreme mass-ratio systems where the 
central black hole is rotating. The motivation for this is clear: Although little is known about the spin distribution 
of astrophysical massive black holes (but see, e.g., [13; 113)' there is no reason to think that massive holes in nature 
are non-rotating. Hence, a useful model of a LISA-relevant inspiral must incorporate a Kerr black hole as a central 
object. Indeed, as this work demonstrates, the spin of the central hole may have a very pronounced effect on the 
value of the self force and hence on the inspiral dynamics. 

Computing the GSF for generic inspiral orbits in Kerr is an extremely challenging task, and this work only represents 
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a first step toward this ultimate goal. The recent advance in calculations of the GSF in Schwarzschild [29f was achieved 
after nearly a decade of development, in which the necessary computational techniques had been devised mainly by 
using the SSF as a simple test bed. In preparing to tackle the Kerr problem, we once again resort here to the 
simplicity of the scalar field toy model. Furthermore, as a primer, we specialize to (geodesic) orbits which are both 
circular and equatorial. This setup already captures much of the complexity of the Kerr problem (and, indeed, 
offers an opportunity to explore some qualitatively new physics), while providing a more manageable environment for 
development. 

Our calculation represents a first application of the standard mode-sum scheme for orbits in Kerr. As such, 
it provides a first test of the regularization parameter values derived in Ref. j30| (we shall review the notion of 
regularization parameters in Sec. Ill below). We opt here to work in the frequency domain, with the obvious advantage 
that we then only need to deal with ordinary differential equations (ODEs). We decompose the scalar field equation 
in a basis of spheroidal harmonics (which are frequency-dependent), and solve the resulting ODEs numerically, with 
suitable boundary conditions. Since the mode-sum scheme requires as input the sp/ierzca^-harmonic modes of the scalar 
field gradient, we then need to re-expand the spheroidal-harmonic solutions into spherical-harmonic components. A 
major technical hurdle intrinsic to this procedure is that the discontinuity of the spherical-harmonic components 
across the particle's orbit hampers the convergence of the frequency series there, due to the Gibbs phenomenon. This 
problem was analyzed in depth in Ref. (sij . and a simple and elegant solution was proposed, which entirely circumvents 
the problem. With this recent development, the frequency-domain approach becomes an attractive option for SSF 
studies, in our view. (We remark that the above Gibbs phenomenon issue does not manifest itself in the case of 
circular orbits considered in our current work.) 

In this work we calculate the dissipative and conservative components of the SSF for a variety of orbital radii and 
black hole spins. Our results for the dissipative component are found to agree well with the numerical results of Gralla 
et al. [Ill (computed from asymptotic fluxes), as well as with the analytic results of Gal'tsov [s^ at large orbital 
radii. As a further important test of our code we verify that the work done by the dissipative component of the SSF 
precisely balances the flux of energy in the scalar waves radiated out to infinity and through the event horizon, as 
extracted from our numerical solutions. For the conservative component our code recovers the results of Diaz-Rivera 
et al. in the Schwarzschild case. This conservative piece is calculated here for the first time for a nonzero Kerr 
spin parameter, revealing several interesting new features. Our main results for the conservative SSF are displayed in 
figure O 

The remainder of this paper is structured as follows. In Sec.|lT]we review the relevant features of circular equatorial 
geodesies of the Kerr geometry, and describe the setup of our problem. In Sec. |lll] we discuss the application of the 
mode-sum scheme for orbits in Kerr, attempted here for the first time. Section IIVI describes our numerical method, 
and in Scc.|V]we provide various validation tests of our code and present our results. Lastly in Sec. IVII we summarize 
our results and consider future work. Throughout this work we use Boyer-Lindquist coordinates (i, r, 6, 0), with metric 
signature ( — h -I--I-) and geometrized units such that the gravitational constant and the speed of light are equal to 
unity. 

II. SETUP AND REVIEW OF PERTURBATION FORMALISM 

A. Orbit and equation of motion 

Consider a pointlike particle of mass ^ and scalar charge 5, set in motion about a Kerr black hole with mass M 
and spin aM . We assume —M < a < M, with negative values of a corresponding to retrograde orbits. We denote the 
particle's worldline (in Boyer-Lindquist coordinates) by x'^{t) and its four- velocity by = dx^/dr, where r is the 
proper time. In this work we neglect the GSF, and consider only the SSF, denoted F^^f{(x q^). Then, the particle's 
motion is governed by 

uf'^piiiu-) = O , (1) 

where the covariant derivative is taken (as elsewhere in this work) with respect to the background Kerr geometry. In 
this work we do not wish to consider the back reaction from the SSF on the particle's motion. Our goal is merely to 
calculate the SSF that would be felt by a particle fixed on a geodesic orbit. We envisage that this SSF information 
could be used to compute the orbital evolution as a second step, but here we do not attempt carry out the evolution 
analysis. For simplicity, we specialize to motion along a geodesic which is both circular [rp(r) = rp =const] and 
equatorial [Op{T) = 7r/2]. Note that, due to the refiective symmetry of the Kerr metric about the equatorial plane, 
an initially equatorial orbit (with 9p = 7r/2 and dOp/dr = at some initial time) would remain so at all times, even 
under the influence of the SSF. 
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Following from the stationarity and axial symmetry of the background Kerr metric, there exist two Killing vectors, 
— dx^/dt and ^^^^ = dx^/d4>. The Kerr metric also admits a Killing tensor Q'^'^ . To each of these there is 

associated a conserved quantity: the specific energy £ = — ^^j-jW^ — ~ut, the specific azimuthal angular momentum 
C = ^'i^^^Un = U0, and the Carter constant Q = Qf^^Uf^Uu. Given initial conditions, these three parameters completely 
specify the orbit of the test particle about the Kerr black hole. 

For our circular and equatorial orbits, one readily finds by solving the geodesic equations (taking 9p — n/2 and 
drp/dr = d'^rp/dr^ = 0) [sl 

l-2v^ + dv^ ^ 1 - 2dv^ + d?v^ , . 

£ = , , C = rnv , , (2) 



where v = \J M /tq and d = a/M . The Carter constant is given explicitly by 

Q^ul + cos^ep[a^{l-£'^)+cBc^epC^] , (3) 
and so it vanishes identically in our case. The angular frequency Vt^ with respect to coordinate time t is given by 



dt g^'i'L - Af (1 + dv"^) ' ^ ' 



where hereafter g^p denotes the Kerr background metric, here evaluated at the circular orbit. Notice our convention 
is that £ and fi^ are always taken positive, with prograde/retrograde orbits distinguished by the sign of a (a > for 
prograde, a < for retrograde). 

Note that in Eq. ([1]) we have kept the mass /i inside the derivative operator. Quinn [l^] (see also Burko et al. j35j ) 
discussed the fact that plausible action principles for the scalar charge in curved spacetime give rise to a dynamically 
varying mass. In general, the evolution of the mass is governed by the SSF component tangent to u": 

'^^-u-F^. (5) 
ar 

In our stationary, circular-orbit setup, however, we must have d^/dr — 0. Therefore u'^Fa — or, more explicitly. 

Ft + n^F^ = 0. (6) 
This trivial relation between Ft and F^l means that in our analysis we need only compute one of these components. 



B. Scalar field equation and multipole decomposition 

We assume that the particle's field $ can be treated as a small perturbation over the fixed Kerr geometry, and that 
it obeys the minimally coupled Klein- Gordon equation 

V„V"$ = -47rr , (7) 

sourced by the particle's scalar charge density T . We model this energy-momentum as a (5-function distribution along 
the particle's worldline, in the form 

T^ql 5\x>^ - x>;{T))[-g{x)]-'l^dT = ^J(r - ro)5(0 - 0p)<5(0 - t^/2) , (8) 

where g = — p^sin^^ is the metric determinant, and where in the second equality we have specialized to Vp = rg 
and 6p — tt/2. The four-velocity component u* is related to the particle's energy and angular momentum through 

Carter discovered fSG*] that the scalar wave equation ([7|) was completely separable in Kerr geometry, with Brill 
et al. giving the explicit separation formula ^37i] . We follow their method and decompose the field into spheroidal 
harmonics and frequency modes in the form 

* = / £ RUnJr)SiJ0; a^)e^"^e-'- du; . (9) 

/=0 m=-/ 
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Here Si^{6; cP') are spheroidal Lengendre functions with (w-dependent) spheroidicity cr^ [we reserve the term spheroidal 
harmonic for the product Sj^{6; (T^)e™"^]. We label the spheroidal Legendre function by Im as we will later introduce 
spherical harmonics which we label by Im. The spheroidal harmonics are orthonormal with normalization given by 

Sirni^-^ <^^y""'^Sp^,{e; a2)e-™'*dr! = 6^,6^^- , (10) 

where the integration is over a 2-sphere t,r =const with area element dQ = sin 0d9d(j)^ and 5„i„2 the standard 
Kronecker delta. 

The source term in Eq. ([7]) is decomposed in a similar manner, writing 



P 



„ OO / 

""T-jT.T. Ti„^J^)SiJ0; a2)e^'"*e-^- do. , (11) 



i=0 m=-; 



where the factor = + cos^ 9 is inserted for later convenience. The periodicity of circular orbits implies that 
the spectrum of the Fourier transform in Eqs. ([9]) and (jlip is given, in our case, by w = nfl^ = uin for integer n. 
Hence for circular equatorial orbits (r^ — tq, Op = 7r/2, (j)p = il^t) T^^^ is given explicitly by 



= |'5,„^W2;a2)<5(r-ro)C, (12) 

where in the second line we have substituted for T from Eq. ([8]). Thus, each m mode contains a single n- harmonic, 
and the spectrum is given by ujn — ujm with 

w,„ = mVL^ . (13) 

Substituting the field and source decompositions into the field equation ^ we subsequently find the radial and 
angular equations to be 

A- (a^%^) + [a'm- - AMrmau^^ + {r^ + a^fu:^ - a^^l^ - A,-„ A)] = -4.Aof^,_,„ (r) , (14) 



dr 

where A = — 2Mr + and Aq = A(ro). The angular equation (|T5t takes the form of the spheroidal Legendre 
equation with spheroidicity — —a^ujf^-^. Its eigenfunctions are the spheroidal Legendre functions S';^(0; — a^cj,^) 
and its eigenvalues are denoted by A^-^. In general there is no closed form for or A^-^^ but they can be calculated 
using the spherical harmonic decomposition method described in appendix [X] When a — the spheroidal harmonics 
S'j-^e*™'^ coincide with their spherical counterparts YJ-^ and their eigenvalues reduce to Aj-^ = l{l + 1). 

As noted by Bardeen et al. (ssj the radial equation (jl4p can be simplified by transforming to a new variable. 



i^imu.Jr)^rR^^^Jr) , (16) 
and introducing the tortoise radial coordinate defined through 

= ? ■ (17) 
dr A ^ ' 

With the above definition the tortoise coordinate is given explicitly in terms of r as 

r.^r^M HNM-') + '» (?^) ' '''' 

where we have specified the constant of integration and r± ~ M ± V M"^ — are the outer and inner roots respectively 
of the equation A = 0. We note that there is an alternative common choice for the tortoise coordinate, namely, 

dn _ r^ + a" , . 

dr " A ' ^ 
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which is useful in that u = t + f* and u = t — r^, are then associated with the "ingoing" and "outgoing" principal null 
congruences of the Kerr background [s^. We shall later refer to in discussing boundary conditions, but for our 
field equation we opt to adopt the coordinate r,, as the coordinate leads to a more complicated radial potential 
[srj . In terms of ipi^^^ (r) and r*, the radial equation (|14p takes the simpler form 



dr2 



47rgAc 



(20) 



where we have substituted for the source from Eq. and where W^^^ is an effective (w-dependent) radial potential 
given by 



ir) 



(r^ + a^)cj„i — am 



A 

l4 



2 2 2{Mr-a^) 
2amijJm + a w H ^ 



(21) 



In the case of circular equatorial orbits, axially-symmetric modes (i.e., ones with m = 0) have vanishing spheroidicity 
/(/ + !). The radial equation (PH)) then admits a simple analytic solution. It is given by 



and A; 



where 



m=0 



airQi{xo)Pl{x) , 
airPi{xo)Qi{x) , 



r < ro , 
r > ro , 



x = I3(r- M) and 



/3: 



M4 



(22) 



(23) 



with xq = a;(ro) and Pj and Q^- being the Legendre polynomials of the first and second kind respectively. The coefficient 
ttj is derived from the jump condition in the derivative of the field at the location of the particle and is given explicitly 
by 



-4^<7(«*/3Ao)-i5,-„(7r/2;0) 



where a prime denotes differentiation with respect to x. 



(24) 



C. Boundary conditions 

Equation (|20|) determines the radial field anywhere outside the black hole once boundary conditions are spec- 
ified on the horizon (r, — > —oo) and at spatial infinity (r* — )■ c»). The boundary conditions follow from physical 
considerations: at the event horizon radiation should be "ingoing" and at spatial infinity radiation should be "out- 
going" (in a sense made precise below). As we approach the boundaries the potential W{r) in the radial equation 
approaches a constant value and the equation becomes that of a simple harmonic oscillator with frequencies 

W^'^ir, ^ oo) = LJr,^ , (25) 

1 /9 . s 2Mr+uJr„ — am 
W^/^ir, -> -oo) = ^ = 7„ . (26) 



Recalling Eq. (|9]) we observe that, at infinity, the t, r-dependence of the /mcj-mode contribution to the full field $ 
will have the asymptotic form exp[— za;m(i ± r»)]/r, where we have converted from to by noting that 

the two coincide (up to an additive constant) at r» — > oo. Choosing the sign such that the exponent is expressed in 
terms of the retarded time coordinate u = t — ensures that any radiation will be purely outgoing at infinity. Hence 
the lower sign applies, and the correct asymptotic boundary condition for the radial field is given by 

^irnJr.^oo)^e-^'^--^' . (27) 

At the horizon the situation is slightly more delicate. The asymptotic radial solutions admit the form ipi„^^ ^ 
exp(±i7mr*) ~ exp[±i(a;m — r7iil+)r»], where we have expressed r* in terms of f* using the asymptotic relation 
r* — > [r+/(2M)]f, + const as r, — oo, and defined 

0+ = — ^ . (28) 
+ 2Mr+ ^ ' 
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(The frequency il+ is the angular velocity fi^ of a stationary observer just outside the event horizon, and might be 
interpreted as the angular velocity of the black hole itself (1^.) In evaluating the fmw-mode contribution to the full 
field $ at the horizon one must now exercise care, and recall that the Boyer-Lindquist coordinate is singular at the 
horizon [s^, and hence the factor exp(im0) in Eq. ^ is singular too. We must instead express the field in terms of 
a regular azimuthal coordinate, and, following [40|, we introduce 

0+ = - n+t. (29) 

In terms of the regular coordinate 0+ we obtain, as r* — — oo, ^{^^^ ~ exp[im(j)+ — i{uJm — mfl+){t =F f*)], where 
=F correspond to ± in the radial solutions ~ exp(±i7m?'*)- For this to represent a purely ingoing radiation the 

lower sign must be selected, so that 3>jj^^ becomes asymptotically a function of only u = t + f , (as well as the regular 
angular coordinates 0+, 6). We thus find that the correct boundary condition at the horizon is given by 

i'imJr. ^ - e-^-- . (30) 

In passing, we remind that frequency modes with ojm < 'm^l^ are superradiant (see, e.g.. Sec. 4.8.2 of [4l|). Since 
in our case Um — mft^, this condition translates to il^ < i7+ [cf. Eq. ([55]) below] and, using Eqs. (|4|) and ([28|) . also 
to ro > rQ{a) where a > and 

/ ^2 X 2/3 

(a) = An^ 



Hence, for prograde circular geodesic orbits with radius greater than rQ'(a), all m-modes of the scalar field are 
superradiant. We will demonstrate this behavior numerically in Sec. IV B] below. 



III. SELF FORCE VIA MODE-SUM REGULARIZATION 

In the standard mode-sum scheme [3, H^l each vectorial component of the SSF is constructed from regularized 
sp/ienca^harmonic contributions, even in the Kerr case. One starts by defining the full force as the field 

Fr'(^) ^ gV„<i>(x) = J2 Pt^'^'i^) ' (32) 

where i^^Q^""-*' denotes the total contribution to Vc($ from its spherical-harmonic Z-mode (summed over m), and x is 
shorthand for x^, an arbitrary field point in the neighbourhood of the particle. Each mode _F'a^""''' is finite at the 
particle's location, although in general the sided limits r — > yield two different values, denoted F^^^"'" respectively. 
The SSF is then obtained using the mode-by-mode regularization formula 

oo C30 

Ff' = E {f^±"^' - - ^ E Pi^'""'^ ' (33) 

where L = I + 1/2 and the regularized contributions Fq'^'^^^ no longer exhibit the ± ambiguity. The (/-independent) 



regularization parameters and Br, were first derived for generic orbits about a Schwarzschild black hole 42| and 
later also for generic orbits in Kerr [30|- In the circular-equatorial orbit case considered here we have At± = Bt = 0, 
and one can show that the mode sum over i^^''^''^^) converges exponentially fast 0. For a — r the regularization 
parameters are generally nonzero and take a rather complicated form; we give these parameters explicitly in appendix 
|B] (specializing to circular equatorial orbits). One usually has f/^'^"^'' cx l~^, so the mode sum in Eq. (1551) converges 
only as ^ l/l. Recall that one can spare the explicit computation of the (f) component F^"'^ by using equation 
Also, from symmetry one obviously has F|°'* = identically. 

In Kerr, as we have seen, the scalar field naturally decomposes into spheroidal harmonic modes and hence in order to 
use the mode-sum scheme in its standard form we must have a preparatory step where the required spherical-harmonic 

modes F^^J^'^-*' are to be constructed out of the spheroidal harmonic modes of the scalar field. To achieve this, we first 
consider the formal expansion of the spheroidal harmonics (with given oj) as a series of spherical harmonics, 

oo 

SiJ0; -a'ioDe^-"^ = E C YUO, 0) , (34) 

i=0 
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FIG. 1: Coupling of spheroidal and spherical modes, illustrated here for a = 0.9M and ro — AM. Shown are the contributions 
from a given [-mode to (left panel) and f^^"'' (right panel), for various spherical harmonic i-modes. (Note that — 
identically for odd values of Z — I.) The width of the I distribution depends mainly on the magnitude of the spheroidicity 
parameter, |(t^| = a?bj^ = a?m?Q.^\ the two cases shown, {l,m) = (44,34) and {l,m) — (44,10), have spheroidicities = 
— 11.821 and —1.022, respectively. The point of this illustration is to note that in practice one only needs to calculate a handful 
more spheroidal I modes than the desired maximum spherical / mode, especially as for smaller a and/or larger ro the coupling 
is weaker than demonstrated above. 



where the couphng coefficients 6|- = are determined as prescribed in appendix [Al [this expansion is similar 

to that apphed by Hughes in Ref."[34| (with a correction noted by Dolan ^43i])]. Note that the spheroidal harmonics and 
the spherical harmonics have the same dependence (i.e., e™"^) and hence only the I modes couple while the m modes 
do not. Using Eq. ^ in combination with Eqs. (fT6|) and (p4| we can then express each of the spherical-harmonic 
l-mode contributions in Eq. p2p (for a = r) in the form 

oo / 

i^(f"")'(x) ^qS/^Y.T. bUiJ^)Yir.{e, 0)e— . (35) 

[=0 m=-i 

The quantities f needed as input for the mode-sum formula (|33p are obtained from the field _F'i^""''' (a;) by taking 
the limits ^ 9p, (f) ^ (f)p and t ^ tp, followed by r — )■ . 

Note in Eq. ((35)) that whilst formally one must sum over all 1 to construct Fq^"""*' , in practice this is not necessary 
as the Z-spectrum (for given /,m) is strongly peaked around I — I; we demonstrate this behavior in figure [TJ The 
bandwidth of / around I increases slowly with increasing spheroidicity |cr^| = a^aj"^, yet even at the largest spheroidicity 
considered in this work (cr^ ^ —126 for a = 0.998M, ro = 2M), we find that only modes within / — 11 < / < / + 11 
carry significant contributions to each of the /-modes fi,^"'''''. 



IV. NUMERICAL IMPLEMENTATION 

For general / and m the radial equation ([20)) has no known analytic solutions so it must be solved numerically. 
To reduce the computation burden one first notes that the individual Im modes of the scalar field are invariant 
under m — >■ — m combined with complex conjugation. Consequently when solving the radial equation we need only 
calculate the modes that have m > as we can recover the negative m modes by taking the complex conjugate of the 
corresponding positive m modes. 



A. Boundary conditions and junction conditions 



The main numerical task is to solve the inhomogeneous radial equation (j20p with the physical boundary conditions 
described by ([27]) and pOp. The form of the inner boundary conditions makes it more natural to adopt as the 
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coordinate for the numerical integration. Our numerical domain extends from r* = r*in <C — M out to r* = r*out ^ M 
(how these boundaries are chosen in practice will be discussed below). We assume that the radial field V'f^ admits an 
asymptotic expansion in l/r at r — >■ cx) and an asymptotic expansion in r — r+ at r ^ r^. . Recalling the leading-order 
behavior of the physical solutions, expressed in Eqs. ^7} and pop . we thus write 



fco 



^ar^ut) = e+— -^crr-';, (36) 

fc=0 

^iJnn) = e-^--'" J2 ^^'(^i" - ^+)' ' (37) 

fc=0 

where Hn = r(r,i,i), Tout = fiTmut) and the truncation parameters fcin,out are chosen such that the boundary conditions 
reach a prescribed accuracy (see discussion below). The expansion coefficients are determined by substituting each of 
the above series into the radial equation. This gives recursion relations for the coefficients c^g'' respectively in terms 
of c^''^^ . These relations are rather unwieldy so we relegate their explicit forms to appendix [Cl 

The homogeneous solutions obtained with the above boundary conditions ([55]) and ([57]) are proportional to the 
yet-to-be-specified constants cj^ and respectively. These constants are determined by imposing suitable matching 
conditions at the location of the particle. The inhomogeneous solution can be written in the form 

^iJ^) = ^,;„(O0('^o -r)+ ^tJr)Q{r - ro) , (38) 

where Q(x) is the Heaviside step function. Substituting this into the radial equation (j20|) and comparing the coefficients 
of the delta function and its derivative we find 



u'Ao 



where a prime denotes d/dr and, recall, Aq = A(ro). The first equation implies that the field is continuous at the 
particle whilst the second describes the nature of the discontinuity in the field's derivative arising from the delta- 
function source. 

In order to determine the correct values of cf^ of eg'*, for which the conditions and PO)) are satisfied, we first 
numerically solve the radial equation (i) starting from the boundary rout with — 1 and integrating inward, and (ii) 
starting from the boundary r-m with Cq'' — 1 and integrating outward. We denote the two corresponding homogeneous 
solutions by tp-' (r) and tp^ (r) respectively, so 

C = ^s°^L ^r™ = ^"^,:;. (41) 

Substituting these relations in Eqs. (15^ and pUj) yields two algebraic equations for cg° and Cg'', whose solutions read 



eh 



^L(^o) 



(42) 
(43) 



Once the coefficients c'^''^^ have been determined, the (unique) physical solution is constructed using Eqs. ([38|) with 

(SD. 



B. Algorithm 



Following is a summary of the numerical procedure we implement for constructing the SSF. We outline the major 
steps and give some details about the numerical method and the choice of numerical parameters. 
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• Fix a black hole spin a and orbit radius rg and calculate the orbital parameters £,C and fi^ [Eqs. ©and ([4])], 
the spherical harmonic decomposition coefficients 6|- and the spheroidal harmonic eigenvalues A^-^ (the latter 

two using the method outlined in appendix for all I and m in the range < I < /max, < m < I. In this 

work we typically take fmax — 55, which is sufficient for calculating all spherical harmonic contributions pj^!^^^^^ 
up to ? ~ 50 in most cases; see below. (The estimation of the contribution to the mode-sum from the remaining 
large-/ tail will be discussed in the next subsection.) 

• For each / mode obtain the axially-symmetric mode of the radial variable, '>Pi m=o^ using the analytic formula 

• (For each m ^ mode) obtain the boundary conditions for the radial variable using Eqs. (1551) and ([57| . setting 
^(x>,eh _ ^ Through experimentation we found it practical to set the inner boundary at r,in = —60 A/. The 
location of the outer boundary required some adjustment depending on the radius of the particle's orbit. In 
practice we took r^,out = 9000M for rg < 30M and steadily moved it outward for increasing rp in order to 
achieve sufficiently fast convergence of the asymptotic series ([551) . The largest value for r+out we used was for 
ro > lOOM where we had to set r^.out = 6.0 x 10''M. We chose fcin,out such that the magnitude of the fcin.out + 1 
term drops below a certain threshold, which we set to 10"^"*. 

• (For each m ^ mode) integrate the homogeneous part of the radial equation (PH]) numerically to obtain (r). 
For this we used the standard Runge-Kutta Prince-Dormand (8,9) method from the GNU Scientific Library 
(GSL) [43|. The GSL Runge-Kutta routine allows one to set a global fractional accuracy target which we took 
here as 10~^^. To test the integrator we used it to solve for a few m = modes and compared with the analytic 
solution (12^ . We made further use of the GSL library to calculate many of the special functions (Legendre 
polynomials, elliptic integrals, Clebsch-Gordan coefficients, etc) that our code requires. 

• Given the numerical solutions -ip^ (For each m ^ mode), proceed to determine the matching coefficients 
^oo,e/i ^.^^ (|42|) and (|43|) . and construct the physical inhomogeneous solutions ip^^ using Eqs. (|4ip and (|38p . 
Record the values of V^^-^ and its (one-sided) r and t derivatives at the radius of the particle. 

• Given V'fm(''o) ^^'^ ^ct±^;"m(^o) for all spheroidal Im modes up to /max, use equation psp to construct the 
spherical-harmonic I modes of the full force at the location of the particle, This procedure allows us to 
obtain all /-modes which do not have significant contributions (through coupling) from the uncalculated modes 
/ > /max- The highest such / mode, denoted /max, is determined by calculating the contributions from the 
/max + 1 spheroidal mode to the various /-modes and identifying the highest value of / for which this 
contribution falls below a given threshold, set here to 10^^^ (fractionally). With /max — 55 we find /max > 44 
for all a, rg within the parameter range considered in this work (lower values of /max for larger \a\ and smaller 
ro, with typical values around /max ~ 50); cf. figure [T] 

• In the final step, calculate the regularized modes Fq'-'^^*-' defined in Eq. (P5)) using the regularization parameters 
given in Appendix [Bl Then sum over / modes as in Eq. (|33p to obtain the desired SSF. Formally, the mode-sum 
formula (I33p requires summation over all / modes from I = to I = oo. In practice, of course, this is neither 
possible nor necessary. For the t component, the mode sum converges exponentially fast, and we typically find 
that the contribution from the modes / > 15 can be safely neglected. For the radial component the situation 
is a little more subtle, as the mode sum converges only as ^ 1// in this case — artificially truncating the series 
at / ~ 50 may potentially result in an error of as much as a few tens of percent in the final SSF. It is therefore 
important to estimate the contribution from the / > /max tail of the mode sum. The method we used for this 
estimation follows that of Barack and Sago [2lJ, and for completeness we review it in the next subsection. 

C. Estimation of the high-/ tail contribution 

We write the total radial component of the SSF as a sum of two pieces, a numerically computed piece, and a large-/ 
tail: 

= ^^<im.x ^^>i„.x ^ (44) 
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where, with fI^'''^^"' as defined in equation (p3 



oo 



_p^<i.„ax = J2 Fl^'^'s) and f;>'»- = ^ F^'-^^e) _ (45) 



1=0 i=i„^x+i 
To evaluate the large-^ tail Ff^^^^^ we extrapolate the last n numerically calculated /-modes using the fitting formula 

AT 



^K.cg) ^ §1 , (46) 

n=l 

where, recall, L = I + 1/2 (how we chose n and N in practice is discussed below). For this fitting we used a standard 
least-squares algorithm from the GSL. Given the coefficients Djn, we then estimate the high-Z contribution using the 
formula 

Af 00 N „^ 

Fl>'--^J2^2n E ^''" = E r9 !"n, ^2n-iamax + 3/2), (47) 

where (x) is the polygamma function of order n defined as 

with T{x) being the standard gamma function. 

Practical use of this estimation method requires some experimentation. For a given N g {3, 4, 5} we considered 
a weighted average of the values obtained for F^°^^ as we vary n from 20 to 35, where the weighting for each term 
is given by the square of the inverse of the fractional difference in the value of F^®'^ as we increase n by one (this 
procedure is meant to bias the average in favour of n values for which F^**'^ depends only weakly on the number of 
fitting modes.) We obtain 3 different average values corresponding to A'' = 3, 4, 5, and use the variance of these values 
to estimate our numerical accuracy (we record as significant figures only those figures that remain fixed as we vary 
N). This error dominates the overall error budget of the SSF, and we hence use it to estimate to overall accuracy of 
our final SSF results. 

It should be noted that the relative contribution from the large / tail is particularly important in the scalar-field 
case (as compared with the gravitational case) . This is because the contribution from the first few I modes turns out 
to be relatively large and opposite in sign with respect to that of the higher modes. In the Schwarzschild case, the 
contributions from the / = 0, 1 modes are both negative and (e.g., for rp = 6M) conspire to nearly cancel out the 
combined contributions from I = 3-6. In the Kerr case this cancellation sometimes involves an even greater number 
of modes (particularly near a, tq values for which the radial SSF vanishes — see below). This behavior is not observed 
in the gravitational case [2l| — at least not for the Lorenz-gauge GSF in Schwarzschild. 



V. CODE VALIDATION AND RESULTS 



A. High-1 behavior 



According to mode-sum theory [l3|, the regularized modes Fr'-'^''^^ in the mode-sum formula ([55)1 should fall off as 
~ 1/P for large I. This behavior relies sensitively on the delicate cancellation of as many as 3 leading terms in the l/l 

expansion of the full modes F^^^"^' (which itself diverges at ~ /), and hence provides an excellent test of validity for our 
numerical results. Indeed, we have been able to confirm a clear ~ 1/P behavior in our numerical data — an example 
is presented in figure [21 Similarly for the time component, we know from theory that the regularized contributions 
pKrcg) (jgpa^y exponentially with I, and again we were able to observe this behavior in our numerical data — see again 
figure [2] for an illustration. The above two tests give us confidence that the high-Z spheroidal contributions (whose 
numerical computation is most demanding) are calculated correctly, and that the spherical-harmonic decomposition 
procedure is implemented properly. These tests also confirm, for the first time, the validity of the regularization 
parameters in the Kerr case (for circular equatorial orbits). 
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FIG. 2: Left panel: the regularized modes Fr'"^"^^' as a function of / for ro — 5M and a — 0.5M. The solid reference line 
is oc l/l^. The regularized modes demonstrate an asymptotic oc l/l^ behavior at large I, as expected from theory (note the 
log-log scale). Right panel: the regularized modes F^^'^'^^^ as a function of I for ro — 5M and a = 0.8M. The solid reference 
line is exponentially decreasing with I. The regularized modes of the t component show a clear exponential decay at large as 
expected from theory (note the semi-log scale). Similar behavior is observed for other values of ro and a. 



B. Energy flux in the scalar waves 



The above validity check only tests the high-/ output of our code. We now discuss a second, more quantitative test, 
which probes primarily the lower-/ portion of the mode sum (and in that sense it is complementary to the first test). 
From global energy conservation we have that the work done by the dissipative piece (here the t component) of the 
SSF must be balanced by the flux of energy carried away in scalar-field radiation. We can use our code to compute 
the flux of energy radiated to infinity and down the black hole, and the result must be consistent with the value of 
the local dissipative SSF. For the t component the mode-sum converges exponentially fast, and it is for this reason 
we argued that the energy-balance test is mostly sensitive to the low-/ portion of the mode-sum. 

We first briefly review the relevant formalism for computing the radiative flux. The stress-energy tensor of the 
scalar field is given by 



ra/3 = ^(<f..<i>,/3-^5o/3<i>'^*M), (49) 



where, as always, gap denotes the Kerr background metric. We wish to consider the flux of scalar-field energy flowing 
to infinity and down the hole. Let E"*" and S~ represent two (timelike) hypersurfaces with r = const 3> M and 
= const <^ —M, respectively; and let dS^ represent a portion of of a small time span dt. The amount of 
scalar-field energy flowing through over time dt is expressed by 

dE± = T T^p^^^^d^t (50) 

(see e.g.. Sec. 4.3.6 of [39!|), where dT,^ represent outward-pointing surface elements on dE^, and the integral is 
performed over the corresponding 2-spheres of constant r, t. The signs are chosen such that the outflow of energy 
through E+ is positive, and so is the inflow of energy through E^ in the Schwarzschild case (recall, however, that 
dE^ can turn negative in the Kerr case, when superradiance is manifest). In coordinate form we have ^^^^ = and 
dE^ = {—g'^^^Y^'^fadOdipdt = S^p^ sin 6 d6d(j)dt, where f;^^-* = —Ap^ sin^ 6 is the determinant of the induced metric on 
E^, and fa = Sai9^^)~^^'^ — SaA~^/^p is an outward-pointing radial vector of a unit length. The (time-independent) 
flux of energy through Ej- is hence given by 

E± = ^^TA<I^Ttrdn. (51) 

From Eq. ([^^ we have Ttr = (47r)^^$_t<i> which, in order to facilitate the angular integration in Eq. (ISlT) . we 
write as Ttr — (47r)~^$_t<i>*^ with an asterisk denoting complex conjugation (this is allowed since $ is a real field). 
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We then substitute the spheroidal-harmonic decomposition 



2, .2 \ ^im<i) ^—iuj^-nt 



(52) 



making the replacement ($f^),t — — *w„i$;-„j. The asymptotic relations 



(53) 



[recall Eqs. p6l) and ([57)) ] also allow us to replace ($^ ),r = — imr^^^^ for r — >■ oo, and ($^ — 2jMr_|-A ^m(ri0 — 
il+)$?^^ for r — > T-i- [where in the last equality we used Eqs. pT|) . ([26)) and ([28)) ]. With these substitutions, the integral 



(54) 
(55) 



in Eq. (|5ip is readily evaluated using the orthonormality relation (|10p , giving 



Ini 



M 



In Table U we display numerical values for the total energy flux, -Etotal = + > ^-s computed using our code 
based on Eqs. ([M)) and ([55)1. For a similar orbital setup, Gralla et al. |32| have previously calculated the total flux 
of scalar-field angular momentum, Ltotai- In the case of circular, equatorial orbits there applies the simple relation 
-Etotai = ^ipLtotai, which allows us a direct comparison with the results obtained in Ref. [33] ■ The data in Table U 
shows good agreement between our fluxes and those of Gralla et ai, with relative differences comparable in magnitude 
to the estimated relative numerical error in the data of Ref. [32] . 

Table U also displays numerical results for the horizon flux, expressed as a fraction of -Etotai- Superradiance 
{E^ < 0) is manifest whenever > fl^. Horizon absorption does not normally exceed ~ 10% even for strong-field 
orbits (as also noted by Hughes [34| in the gravitational case), but prograde orbits around a fast rotating hole can 
display extreme superradiance behavior [nearly 25% negative absorption in the example of (a, tq) = {0.998M,2M)]. 
The graph in figure ([3]) displays some more horizon absorption data. 
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FIG. 3: The horizon flux of scalar-field energy, E-, as a percentage of the total fiux for different orbital radii ro and spin 
parameters a. The curves are interpolations based on the numerical data points shown. Superradiance behavior {E- < 0) is 
manifest whenever the horizon's angular velocity f2+ is greater than that of the particle. 
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a/M 


ro/M 






E-fEtotal 1 


— iitotal/ -C/total 


1 — -Etotal/f 


0.998 


2 


4.3975979e- 


-3 


-0.2486 


7.06e-7 


-4.7e-10 




4 


6.65618888e 


-4 


-0.1168 


2.12e-7 


-1.6e-10 




6 


1.69712483e 


~4 


-0.0692 


1.12e-6 


-9.2e-ll 




8 


6.04494314e 


-5 


-0.0464 


-2.12e-6 


-4.6e-ll 




10 


2.64845608e 


-5 


-0.0337 


6.65e-8 


-3.7e-ll 




20 


1.87388789e 


-6 


-0.0120 




-2.8e-12 




40 


1.23796212e 


-7 


-0.0041 




7.7e-ll 


0.5 


6 


2.02918608e 


~4 


-0.0248 


-5.19e-7 


-8.9e-ll 




8 


6.76202950e 


-5 


-0.0196 


-1.76e-6 


-6.8e-ll 




10 


2.86637838e 


-5 


-0.0151 


7.33e-7 


-3.3e-ll 




20 


1.92605066e 


-6 


-0.0058 




-l.Oe-12 




40 


1.24998716e 


-7 


-0.0021 




-3.5e-ll 


0.0 


6 


2.55199967e 


-4 


0.0308 




-9.2e-ll 




8 


7.72547978e 


-5 


0.0114 


1.98e-6 


-6.2e-ll 




10 


3.13766525e 


-5 


0.0054 


1.28e-7 


-4.1e-ll 




20 


1.98366995e 


-6 


0.0006 




-4.6e-12 




40 


1.26226716e 


-7 


0.0001 




4.2e-ll 


-0.5 


8 


9.02315446e 


-5 


0.0468 


-5.01e-7 


-4.9e-ll 




10 


3.47579647e 


-5 


0.0284 


5.20e-6 


-4.6e-ll 




20 


2.04718763e 


-6 


0.0073 




3.4e-12 




40 


1.27600490e 


-7 


0.0022 




5.2e-ll 


-0.998 


9 


6.22560292e 


-5 


0.0644 


-7.86e-7 


-5.0e-ll 




10 


3.88839360e 


-5 


0.0519 


-1.56e-6 


-4.2e-ll 




20 


2.11643277e 


-6 


0.0142 




2.2e-ll 




40 


1.28992555e 


-7 


0.0044 




-6.1e-ll 



TABLE I: Scalar-field energy flux for various values of the spin parameter a and orbital radius ro. The 3rd column displays 
the total fiux of energy radiated to infinity and down the black hole, as extracted from our numerical solutions. The 4th 
column presents the fraction of the total power absorbed by the black hole, with negative values indicating superradiance. 
The 5th column compares our fluxes to those obtained by Gralla, Friedman and Wiseman (GFW) [33], showing a good 
agreement. (GFW provide results for the radiated angular momentum, which we convert here to radiated energy using the 
relation -Btotai = fi^itotai; their results are given with 6 significant figures.) In the last column we test our SSF results (for 
the dissipative component) against the balance relation (|57p as discussed in Sec. IV Cl £{< 0) is the rate at which the particle's 
scalar energy is dissipated, as computed from the local SSF using Eq. (I56|l . In this Table (and all subsequent Tables) we use 
an exponential notation whereby (e.g.) 'e— 3' stands for xlO"''. All decimal places presented are significant. 



C. Dissipative component of the SSF 

In the case of circular, equatorial orbits the entire information about the dissipative effect of the SSF in contained 
in the two components Ft and F^. Specifically, we obtain from Eq. ([1]) 

fi£ = -{uY^Ft, lit = {u*)-^F^, (56) 

where, as elsewhere in this work, an overdot denotes d/dt. The relation ^ implies that in practice we need only 
calculate one of the two components Ft and F^ — here we choose to calculate the former. Sample numerical data for 
Ft are presented in table HIl 

In our stationary setting, the rate at which the particle is loosing scalar energy, given by must equal the rate 
at which energy flows to infinity and down the black hole, given by i^totai. Using Eq. (|56p we may express this energy 
balance relation directly in terms of the SSF: 

Ft = -fiu^S = /iu*£'totai • (57) 

As discussed above, this allows us to test our computation of Ft (primarily the low-/ portion of the mode-sum) 
by verifying that our numerical results satisfy Eq. (j57p . As the data presented in right-most column of Table |T] 
demonstrate, we indeed find a very good agreement. 
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{M^/q^)Ft 


ro/M 


a = -0.9M 




a = -0.7M 






a = -0.5M 




a = 






a = O.SAf 






a = 0.7M 






a = 0.9M 




4 


_ 




_ 






_ 




_ 






_ 




1 


359218156- 


3 


1 


142048206- 


3 


5 


- 




- 






- 




- 




6 


07684087e- 


4 


5 


357685616- 


4 


4 


796349856- 


4 


6 


- 




- 






- 




3.60907254e- 


4 


2 


783947986- 


4 


2 


551610136- 


4 


2 


357338536- 


4 


7 


- 




- 






- 




1. 767320196- 


4 


1 


46366447e- 


4 


1 


371037036- 


4 


1 


290467476- 


4 


8 


- 




- 




1 


15781360e- 


4 


9.77204485e- 


5 


8 


44876316e- 


5 


8 


024073936- 


5 


7 


645191606- 


5 


10 


4.60475173e- 


5 


4.38590519e- 


5 


4 


18429073e- 


5 


3.75022727e- 


5 


3 


40410532e- 


5 


3 


286111976- 


5 


3 


177601686- 


5 


14 


1.03173965e- 


5 


1.00539090e- 


5 


9 


80387438e- 


6 


9.23672660e- 


6 


8 


74728207e- 


6 


8 


570772246- 


6 


8 


403735786- 


6 


20 


2.28457108e- 


6 


2.25311511e- 


6 


2 


22274047e- 


6 


2.15159216e- 


6 


2 


08709237e- 


6 


2 


063009026- 


6 


2 


039805746- 


6 


30 


4.30761267e- 


7 


4.27729235e- 


7 


4 


24767592e- 


7 


4.17678576e- 


7 


4 


11035602e- 


7 


4 


084969126- 


7 


4 


060210076- 


7 


50 


5.43419839e- 


8 


5.41729302e- 


8 


5 


40064364e- 


8 


5.36016621e- 


8 


5 


32132722e- 


8 


5 


306236476- 


8 


5 


291388076- 


8 


70 


1.40256823e- 


8 


1.39999178e- 


8 


1 


39744575e- 


8 


1.39121644e- 


8 


1 


38518165e- 


8 


1 


382821036- 


8 


1 


380489826- 


8 


100 


3.35072295e- 


9 


3.34717963e- 


9 


3 


34366914e- 


9 


3.33503895e- 


9 


3 


32661812e- 


9 


3 


323307556- 


9 


3 


320029176- 


9 



TABLE II: Sample numerical results for the t component of the SSF. Entries left empty correspond to orbits below the 
inner- most stable circular orbit (ISCO). All figures presented are significant. 




FIG. 4: Time component of the SSF; comparison with Gal'tsov's slow-motion formula. Plotted is the relative difference between 
our "full" SSF Ft and Gal'tsov's weak- field/slow motion analytic approximation (|58l) as a function of orbital radius tq. Solid 
lines are interpolations of the data points shown. We show results for a — ±0.998M; similar agreement between Ft and f{^^^ 
at large ro is manifest for other values of a too. 



It is also interesting to test our results against the vifeak- field/slow-motion analytic formula derived by Gal'tsov '33^, 

^Gartsov ^ 1^2^^ (^^2^3 ^ - , (58) 

which is valid for tq 3> M. Here the first term corresponds to the radiation heading out to the infinity and the second 
to the radiation absorbed by the black hole. In figure|4]we plot the relative difference between the "full" SSF computed 
here and ir-Gai'tsov function of tq for a couple of a values (we choose the two extreme cases a — ±0.998M). Our 
results seem to obey Gal'tsov's formula for large orbital radii, as expected. 

Lastly, we note that our value of Ft for (a, tq) = (0, 6Af) (see Table HI)) coincides through all 9 significant figures 
with the value computed by Haas and Poisson in ReL |45|. 
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D. Conservative component of the SSF 

In our orbital setting, the conservative effect of the SSF is entirely accounted for by its radial component, Fr- The 
computation of this component is more involved, as in this case the mode-sum requires regularization, and (relatedly) 
the mode-sum series exhibits slow convergence. While results for the dissipative SSF in Kerr (obtained indirectly 
from the asymptotic fluxes) already exist in the literature, our results for Fr are new. 



{M^/q^)Fr 


ro/M 


a = -0.9M 


a = -0.7M 


a = -0.5M 


a = 


a = 0.5M 




a = 0.7M 




a = 0.9M 


4 
















-5.24194e- 


4 


-9.5941e- 


4 


5 












-4.160235e- 


-5 


-2.044174e- 


-4 


-3.63448e- 


-4 


6 










1.677283e-4 


-2.421685e- 


-5 


-9.528095e- 


-5 


"1.645525e 


-4 


7 










7.850679e-5 


-1.467677e- 


-5 


-4.980678e- 


-5 


-8.410331e 


-5 


8 








9.642777e-5 


4.082502e-5 


-9.21907e- 


6 


-2.829488e- 


-5 


-4.696081e 


-5 


10 


4.939995e-5 


4.100712e 


-5 


3.28942e-5 


1.378448e-5 


-4.03517e- 


6 


-1.091819e- 


-5 


-1.768232e 


~5 


14 


9.968208e-6 8.303689e 


-6 


6.67043e-6 


2.720083e-6 


-1.07573e- 


6 


-2.561183e- 


-6 


-4.02935e- 


-6 


20 


1.878548e-6 


1.565128e 


-6 


1.2550019e-6 


4.93790e-7 


-2.50260e- 


7 


-5.43942e- 


7 


-8.35474e- 


-7 


30 


2.873310e-7 2.389538e 


-7 


1.90843e-7 


7.1719e-8 


-4.595209e- 


-8 


-9.26682e- 


8 


-1.391883e 


-7 


50 


2.74358e-8 


2.272902e 


-8 


1.803392e-8 


6.3467e-9 


-5.27419e- 


9 


-9.90589e- 


9 


-1.452810e 


-8 


70 


5.87543e-9 


4.8525e- 


9 


3.8312e-9 


1.2845e-9 


-1.25352e- 


9 


-2.26649e- 


9 


-3.27820e- 


-9 


100 


1.1508e-9 


9.4715e- 


10 


7.4364e-10 


2.356e-10 


-2.7134e-10 


-4.7388e-10 


-6.7625e- 


10 



TABLE III: Sample numerical results for the r component of the SSF. Entries left empty correspond to orbits below the ISCO. 
All figures presented are significant. The numerical accuracy is lower compared to that of Ft as a result of (i) the regularization 
procedure involved in obtaining Fr, and (ii) the slow decay of the large-/ tail in the case of Fr (compared with the exponential 
decay of the tail for Ft). 



Table Hill presents Fr data obtained for a range of a and tq values. Our results for Schwarzschild (a = 0) agree with 
those of Diaz-Rivera et al. [l3| through all significant figures. The most striking feature of our results is that — unlike 
in the Schwarzschild case where the radial SSF is always repulsive (outward pointing) — here we find that for certain 
prograde orbits Fr becomes attractive (inward pointing). This behavior is better illustrated in figure [5l where we 
present a contour plot of Fr across the parameter space of a,rQ. This plot is based on the data shown in table Hill 
as well as many other intermediate data points. A few fixed-rg and fixed-a cross-sections of the contour plot are 
presented in figure [S] for clarity. 

We observe the following: (i) For retrograde orbits (a < 0) the radial SSF is always repulsive, as in the Schwarzschild 
case, (ii) For prograde orbits (a > 0) there exists an a-dependent radius Tc at which the radial SSF vanishes; it is 
repulsive for vq < and attractive for rg > Tc- (iii) The critical radius r,. decreases monotonically with increasing 
a. (iv) The critical orbit coincides with the ISCO for a ~ 0.461M; hence, all stable circular geodesies experience an 
attractive radial SSF when a > 0.461Af. It is interesting to note that Burko Tg*] observed a similar change of sign in 
the radial SSF when studying accelerated (non-geodesic) circular orbits in Schwarzschild geometry. 

To gain some intuition about the above behavior of the radial SSF, it is instructive to analyze our results in the 
context of post-Newtonian (FN) theory. In the Schwarzschild case, a weak-field expression for the radial SSF was 
worked out to high FN order by Hikida et al. in Ref . ^] . Only the leading 3PN and 4PN terms are given explicitly 
in that work. They read [49| 



P3+P3°^ln(ro/Af) 



P4+Pi^Hro/M) 



(59) 
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50 60 

ro/M 



FIG. 5: The radial component of the SSF, multiphed by rp for convenience, across the a, ro parameter space. Contour hnes 
are hnes of fixed r^Fr, with labels giving the value of {M / q)^ {rg / M)^ Fr ■ The near-vertical thick line indicates the location of 
the ISCO, while the near-horizontal thick line marks the curve ro = rc(a) along which the radial SSF vanishes. The two lines 
intersect at a ~ 0.461M; for a > 0.461A/ all stable circular geodesies experience an attractive radial SSF. 




FIG. 6: Left panel: Radial component of the SSF as a function of a for various fixed values of the orbital radius ro. Right 
panel: Radial component of the SSF as a function of ro for various fixed values of the spin parameter a. In both panels dots 
represent numerical data points, and solid lines are interpolations. 



where the coefficient are given by 



3 64 9 



P3 



14 66, 29 , 604 

P. = -y^"Tl-2 + ^.^ + -. 1.85852..., 

pr = I , (60) 

with 7 = 0.577215 . . . being the Euler number. Note the leading 3PN term is dominated by a "logarithmic running" 
term. Using Eq. (j59p as an ansatz for a = 0, we performed a two-dimensional fit of a large-rg subset of our numerical 
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data to a model of the form Fr = Fr + aCx power series in M/rg. We find, at leading order, 

Fr{r » M) = +pT^ f-V , (61) 

with 

pI° ~ -1.00091. (62) 

Our numerical accuracy was not sufficient to distinguish between different PN models (including possible logarithmic 
terms) at higher PN orders, so we do not present here fit results beyond the leading 3PN spin term. This leading 
term has the interpretation of a spin-orbit coupling ( "a • i" ) . We are not aware of any explicit analytic calculation 
of this term in the PN literature. (It might be possible to extract the 3PN spin-orbit term from the formal results 
of Ref. [13], which, however, we have not attempted here.) Our numerical fit suggests that the coefficient p|° of the 
leading 3PN spin-orbit term is simply —1. 

In figure [7] we plot some of our Fr numerical data points against the analytic PN model (|61|) . A good agreement 
is manifest down to radii as small as tq = lOM where the difference between our fitted PN formula (pT|) and our 
numerical results is in all cases no more than 8%. At ro — 20 M this difference is never greater than 3%. 



0.0001 




ro/M 

FIG. 7: Comparison of numerical data for Fr (dots) with the PN fit model (|6ip (solid lines). For prograde orbits with 
a < 0.461Af the radial SSF changes sign at ro — rc{a); cf. figures [5] and [6] 

We note that £ ~ Tq' for large ro [recall Eq. ([2|)], and hence the leading spin term in Eq. ([6T1) dominates the overall 
behavior of Fr at sufficiently large rp, falling off as ^ 'Tq'*"^. At intermediate values of ro, this term, which is negative 
for a > 0, competes with the leading "Schwarzschild" term, which falls of as ^ r^^ Inrg and is positive. This, we now 
observe, gives rise to the change-of-sign observed for Fr in our numerical data. 

VI. CONCLUDING REMARKS AND FUTURE WORK 

In this work we presented a first calculation of the SSF experienced by a particle orbiting a Kerr black hole, 
specializing to circular and equatorial geodesic orbits. This represented a first application of the mode-sum method 
in Kerr, and as_a by-product we confirmed the analytic values of the regularization parameters Aa, Ba and Cq,, as 
calculated in [30j, for the above class of orbits. Our numerical calculation relied on a standard frequency-domain 
decomposition of the scalar field equation in terms of spheroidal harmonics; the spherical-harmonic contributions 
required within the regularization procedure were obtained by projecting the spheroidal-harmonic contributions onto 
a basis of spherical harmonics. 
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We tested the performance of our code in various ways. The contribution to the SSF from the high-Z modes was 
found to possess the expected behavior, falhng off exponentially for the time component and as ~ for the radial 
component. We confirmed numerically that the work done by the time component of the SSF precisely balances the 
energy in scalar waves radiated out to infinity and down through the event horizon. The energy flux calculated from 
our code also agreed closely with the previous numerical results by Gralla et al. [s^l as well as with Galt'sov's analytic 
formula [s^ in the large-ro regime. The radial, conservative component of the SSF was calculated here for the first 
time. Our code produces good agreement with the previous results of Diaz-Rivera et al. flT'l in the Schwarzschild 
case. For non-zero spin, we observed a qualitatively new behavior: The SSF on prograde orbits with radius larger 
than a certain a-dependent radius Tc turns from repulsive (as in the Schwarzschild case) to attractive. While we 
have no genuine physical intuition to explain the direction of the radial SSF (not even in the Schwarzschild case), we 
observed, at a formal level, that the above change-of-sign may be attributed to a competition between a repulsive 
"Schwarzschild" term and an attractive spin-orbit coupling term. 

This observation came from fitting our numerical SSF data to an analytic FN model at large tq. We thus derived 
a numerical approximation for the leading-order, 3PN spin-correction term. It would be interesting to test our result 
against an analytic PN computation of the radial SSF, once the FN result becomes available. To further make contact 
with PN theory it would be necessary to extract higher-order terms in the PN series, and for this it may be necessary 
to improve the accuracy of our code at large orbital radii. The main limiting factor, and by far the dominant source of 
error in our calculation, is the large contribution to the SSF from the long uncomputed tail of the /-mode series. The 
relative contribution of this tail increases with tq; in our analysis the uncomputed tail contribution for ~ lOOM is 
more than twice that of the computed modes! The problem can be mitigated in future work by pushing our numerical 
calculation to higher Z, or — better still — by obtaining analytic expressions for some of the higher-order terms in the 
l/l mode-sum, thereby accelerating the convergence of the mode-sum. (This latter technique was applied successfully 
by Detweiler et al. in the Schwarzschild case [48|.) 

As mentioned in the introduction, in general a frequency-domain application of the mode-sum method is made 
difhcult by the bad convergence of the frequency mode sum along the particle's orbit ("Gibbs phenomenon"). The 
problem is unnoticed for circular orbits, since in this case the scalar field is a smooth function of time along the orbit. 
However, the issue will need to be addressed in contemplating the extension of our code to more generic orbits. The 
recently introduced method of "extended homogeneous solutions" 31] proposes a simple method to overcome the 
above difficulty and we envisage incorporating this method in a future extended version of our code. We have already 
started work to generalize the code to eccentric orbits (which, as a first step, we keep equatorial). 

Extension to the gravitational problem is more challenging. The main obstacle is the lack of a formal framework 
for analyzing Lorenz-gauge metric perturbations in the frequency-domain in Kerr. A potential avenue of approach 
would be to work with coupled tensorial spherical-harmonics, although this may pose a significant technical challenge. 
Another possibility would be to develop a suitable tensorial spheroidal-harmonic basis for decomposition in Kerr, akin 
to the tensorial spherical harmonics that can be used in the Schwarzschild case. 
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Appendix A: Spheroidal harmonics and their expansion in spherical harmonics 

The spheroidal harmonics ^^-^^(f?; (J^)e""'^ satisfy the differential equation 

%n(^;^')e™* = 0, (Al) 

where the constant parameter is the spheroidicity. The functions Si^^{9;a'^)e"^'^ are called oblate or prolate 
spheroidal harmonics, depending on whether cr^ is negative or positive, respectively. A useful and efficient method 
for calculating the spheroidal harmonics is via decomposition in spherical harmonics. This method is doubly useful 
in our case, as it automatically generates the spherical-harmonic data required as input for the mode-sum formula. 
The expansion of a given spheroidal harmonic as a series of spherical harmonics, for given m, takes the form 

oo 

Si^^{0;a')e"^^ = E b\j<j')Yi^{e, c^) , (A2) 
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where Imin = \^\- In order to calculate the coefRcients b\ we substitute this expansion into equation (jAip . Noting 



that the satisfy (jAip when cr = with A;,„ = l{l + 1), we get 

oo oo 

^ b\Ja' cos^ + + /)]1^,„ = ^ 6;-^ll„, . (A3) 

Next we multiply the above expression by Y~* and integrate over the sphere. The resulting inner products are given 
by 

'Vfin^dn = , (A4) 



1 2 /2/ + 1 



'imdn = -5^^ + -^{l,2,m,0\l,m){l,2,0,0\l,0) = ki^ . (A5) 

Here the numbers {ji , :/2 , , Ijm) are standard Clebsch-Gordan coefRcients, the form of which implies that k'- ^ 
only for ^ G {/ — 2, / — 1, Z, ^ + 1, ^ + 2}. Consequently, Eq. (IA3p reduces to the recursion relation 

a'd-^l/r^ + a'kl-'b'r^ + [a'kl + l{l + 1)]M + a^l^+%+^ + a^}^+%\+^ - A. 6^ (A6) 

/m im Im Im /m ^ " Im Im Im Im Im im ^ ' 

for the expansion coefficients (with given l,m). This can be put in a matrix form, Kh ~ Ab (keeping the indices 
I, m implicit), where is a known band-diagonal matrix (made up of the known and fe[^) and b = {b'j^, ^^{^ i ■ ■ ■)■ 
This is a standard eigenvalue problem for the eigenvectors b and eigenvalues A (for each m), and the band-diagonality 



spheroidal-harmonic eigenvalues A,- , which we adopt in this work, follows closely that of Hughes in [34 



of K makes is readily amenable to numerical treatment. This method of obtaining the expansion coefRcients h~ and 



Appendix B: Regularization parameters in Kerr geometry 

The regularization parameters for the SSF in a generic orbit about a Kerr black hole were calculated by Barack 
and Ori and given in their Ref. [13 (see for a detailed derivation). For circular, equatorial orbits they reduce to 

- = , (Bl) 

and (in Boyer-Lindquist coordinates) 

= T'Z^A-i/2(g^^ + /:2)-i/2 (32) 

Af = A± = A± = , (B3) 

where the metric function g^^ is evaluated on the equatorial orbit. The expression for is more complicated. It 
can be written in the form 

- q^P^abcdl'''"' , (B4) 

where hereafter Roman indices run over the two Boyer-Lindquist angular coordinates 0, 4> only. The coefRcients Pp,abcd 
are given by 

Pt^abcd = {^TT)-^[iPt.dPabc - (2P^afc + Pab^.)Pcd] , (B5) 

where 

PaP = 9al3 + UaUp , (B6) 

PaPy = (uxUjT'^p + gap,'y/2) , (B7) 



with the Kerr connections and metric functions gap all evaluated on the equatorial orbit. The quantities I"^'"^"- 
are 

jabcd^ / G(7)-^/2(sin7)^(cos7)^-^d7 , (B8) 
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where 

G(7) = sin^ 7 + 2Pe^ sin 7 cos 7 + Pee cos^ 7 , (B9) 
and N = N(abcd) is the number of times the index <j) occurs in the combination (a, 6, c, d), namely 

iV = <5^ + <5^ + 55 + 5^. (BIO) 

The quantities can be written expUcitly in terms of complete elliptic integrals 0, [s^]- In the case of a circular, 
equatorial orbit these expressions become 

jabc, ^ 2(i-^)/WxH + /W£:H 

24P;/V(l_^„)2 

where K{w) = Jq^'^{^ — w siii^ x)~^^'^ dx and E{w) = /J^^^(l — w sin"^ x)^^'^ dx are complete elliptic integrals of the 
first and second kind respectively, and 

u,^l-^. (B12) 

The coefficients l]^'' and I^'' are given by 

iP = 16u;2(2 - 3w) , 1^ = 64^2 (2w - 1) , 

4P=4^^=0, lP=32wHw-l), 

if = 32wHw' -3w + 2), if = if = , ^ 

if = -16^2(^2 + w - 2) , if = -64^2(^3 -w^\) . 



Appendix C: Boundary conditions for the radial scalar- field equation 

In order to derive recurrence relations for the asymptotic expansion coefficients Cj°° and c^'* in Eqs. ([55]) and p7l) . 
we substitute these equations into the homogeneous part of the radial equation ([20]) . By comparing the coefficients of 
(at infinity) or (r — r^)^ (at the event horizon) we obtain 5-term recurrence relations for each of and c|^^q. 

Setting c^g^ = and c^q'' = 1 determines all coefficients c^g'' in a recursive fashion. 
Explicitly, the above recurrence relations are given by 

5 5 

where the various coefficients f°° and /,f read 

= -2A:w™i, 

= fc2 _ ^.^^ ^ cj,„(a2^^ „ 4iM) + k{4iMujm - 1) , 

/2°° = 2[m2(2 - k)ojrn + M {a^Lol, - 2amw„ - 2fc2 + 5fc - 3 + Aj^)] , 

/3°° = 4(fc - 2)2m2 - a'i^^^ - 2fc2 + 8fc - 8 - m2) , 

= -2a2Af (2/c2 - 11/c + 15) , 

= a* (A:2 - 7fc + 12) , (C2) 
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f^'' = -a* [fc^ + fc(-3 - ^ir+LJm) - r^w^ + 2] + 2a^mr+{r+Um + 2ik) 

-c?r+ [2M (6fc2 + A;(-12 - 6ir+a;„) - r^w^ + 3) + r+ (-12^^ + 2fc(9 + 8ir+a;„) + r^w^ + A^-^ - 2)] 

+2oTOr^ [r+a;™(r+ - 2M) - 2ik{SM - 2r+)] 

+4 [4 - 9fc + 1) + 2Mr+ (-20A;2 + lOifcr+w^ + 24A; + A,-^ - l) 
+4 (l5/c2 - Uikr+LOm - 15/c - A;„)] , 
fi^ = 2 {a^LOm{ik + r+u!„i — i) — ia'^m{k — 2ir+uJm — 1) + a^M [2fc^ — /c(9 + 6zr+Wrn,) — Sr^to^ + Gir+cOm + lO] 
Wr+ [-4k^ + 3A;(5 + 4ir+w„) + 2r+w^ - 12zr+w™ + A,-^ - 13] 
+2amr+[3M{ik + r+cOm — i) + r+{—3ik — 2r+LJm + 3i)] 

+r+ (-8fc^ + 30k - 26) + Mr^ (20/;;^ - 20i/;;r+a;„ - 66A; + 20ir+a;„ - 3A,-^ + 49) 

+rl (-lOfc^ + 15ikr+LJm + 30k - ISir+w^ + 2A;-^ - 20)} , 
/I'* = a^u;^-2a3mw,„ 

+0^ [-2fe2 + 4k{-iMujrn + 4ir+LUrn + 3) - 6Mr+w^ + 8iMw„ + 6r^w^ - 32ir+w„ + A;-^ - 18] 

+Aiam[M{k - 3ir+ujm - 2) + r+{-2k + 3zr+Wm + 4)] 

+2Mr+ (lOfc^ - 20ikr+LJm - 54A; + 40ir+a;„ - 3Aj^ + 7l) 

(-15fc2 + AOikr+uJm + 75k - SOir+Urn + 6Aj^ - 90) - 4{k - 3fM'^ , 
= -2a^u!m{-2ik + Mujm - 2r+uJm + 6i) + 2am{-ik + 2Mu)m - ^r+ujm + 3i) 

+M [4fc2 - fc(30 + 20ir+LOm) + 60ir+a;„ - 2A;^ + 56] 

+2r+ [-3fc2 + 3fc(7 + 5ir+uj^) - 45ir+cj„ + 2A;^ - 36] , 
ft = a^uol, - 2amu}m - k"^ + k{-4iMuJm + 12ir+w„ + 9) + l&iMuJm - 48ir+w„ + A^-^ - 20 , 
ff = 2i{k - 5)uJm ■ (C3) 
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